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43 ■ Abstract 

By employing an analytically solvable model including the Duschinsky rotation effect, we in- 



vestigated the applicability of the commonly used Born-Oppenheimer (BO) approximation for 
separating the proton and proton donor-acceptor motions in theories of proton coupled electron 
transfer (PCET) reactions. Comparison with theories based on the BO approximation shows that, 
the BO approximation for the proton coordinate is generally valid while some further approxima- 
tions may become inaccurate in certain range of parameters. We have also investigated the effect 
of vibrationally coherent tunneling in the case of small reorganization energy, and shown that it 
plays an important role on the rate constant and kinetic isotope effect. 
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I. INTRODUCTION 



Proton coupled electron transfer (PCET) reactions are important in many chemical, 
biological and electro-chemical processes. l|-|5| We consider in this paper the case of con- 
certed PCET reactions, where the proton and electron are transferred simultaneously with 
no stable intermediates. Theories of PCET reactions were developed over the past two 
decades. Cukier has proposed a theory for concerted PCET based on nonadiabatic tran- 
sitions between multiple vibronic states, p-[8| which treats the proton coordinate as a high 
frequency intramolecular vibrational mode, and leads to rate constant expressions similar 



to the Bixon-Jortner model [9] for the electron transfer (ET) reactions. |lOl. Ill| Hammes- 
schiffer and coworkers have generalized this theory to include the environmental collective 
coordinate that couples to the proton coordinate. 12Ml4J| 

Recent works on the PCET theory have focused on the importance of the proton donor- 
acceptor motion. The main concept is that fluctuations of the donor-acceptor separation 
Rda strongly affect the vibrational wave function overlap for the proton coordinate, and thus 
the effective transfer integral between vibronic states. Such effect has been first investigated 
in the case of vibrationally nonadiabatic proton transfer reactions, and is found to play an 
important role and leads to much smaller kinetic isotope effect (KIE) than simple estimation 
using the overlap of vibrational wave functions at the equilibrium proton donor-acceptor 
distance. 15Ml8| The effect of the donor-acceptor motion has also been widely discussed in 



recent PCET theories. 



19 



25[ In these theoretical treatments, the proton degree of freedom 
is first quantized, and a Born-Oppenheimer approximation is applied to separate the proton 
motion and the slower degrees of freedom that couple to it. The PCET problem is then 
reduced to transitions between a group of vibronic states. Within the above theoretical 
framework, the effect of donor-acceptor motion can be treated either statically by doing a 



thermal average over the Rda distribution 
a time correlation function formalism. 20-22 



23 



25] 



24| or can be treated dynamically using 



It is interesting to investigate the applicability of the BO approximation to the proton 
motion for several reasons. First, although the masses of proton and deuterium are much 
less than those of the heavy atoms, the difference is not as dramatic as the mass ratio 
between the electron and nuclei. When the hydrogen bond is stiff between the donor and 
acceptor, high frequency donor-acceptor motion raise the question whether a separation of 



time scales is still appropriate. Second, even after the BO approximation is applied to sep- 
arate the proton motion and the proton donor-acceptor motion, additional approximations 
such as static donor-acceptor motion, 



19 



23 



24| or analytical approximations on the vibra- 



tional wave function overlap are often employed. 20|, |22|, |25| It is also desirable to test these 



approximations quantitatively using a solvable model. 

In this paper, we provide tests of the above mentioned BO approximation for the proton 
motion using a model Hamiltonian for concerted PCET reactions. We show that, under 
certain well defined approximations, the effect of proton donor-acceptor motion can be in- 
corporated in a model Hamiltonian base on the Duschinsky rotation effect (DRE). [26] The 
DRE describes the mixing of the normal modes between the donor and acceptor poten- 

I L 

tial energy surfaces. [261] which is different from the linear displaced harmonic modes used 



in the conventional spin-boson model. 27 



During the past years, the effect of DRE 



have been widely discussed in areas such as electronic spectroscopy [29N32J nonadiabatic 



relaxations. 



29 



33 



34j and electron transfer reactions. 35 



39J However, to the best of our 



knowledge, a model based on the DRE has not been applied in studies of PCET previously. 
When the electronic coupling is small, the DRE model allows us to calculate the PCET 
reaction rates analytically using the Fermi's golden rule. This analytically solvable model 
is then applied to assess the BO approximation for the proton motion, as well as the static 
donor-acceptor motion approximation and the analytical approximations to the vibrational 
wave function overlaps. The remainder of this paper is arranged as follows. The model 
Hamiltonian and theories to calculate PCET rate constants are presented in Sec. [Til The 
results are presented in Sec. IHIl where PCET rates and the KIEs are calculated using the 
FGR, and compared with results obtained from approximate theories. We have also studiec 



the rate constant and KIE in the case of small electronic reorganization energy happen, 25 



40 



42| | with the focus on the possible role of vibrationally coherent tunneling. 43|, |44J The 



main conclusions are summarized in Sec. IIVI 



II. THEORY 

A. Model Hamiltonian 

We consider the case of concerted proton and electron transfer with a small electronic 
coupling (nonadiabatic ET), where a two electronic surfaces description is sufficient. |l4. 



22| and the Fermi's golden rule (FGR) that treats the electronic coupling to first order 



perturbation is valid. 



To this end, we apply an extended spin-boson model to describe 



the PCET reactions. The total Hamiltonian can be written as 

H = H S + H B + H BS , (1) 

where, the system Hamiltonian H$ includes terms related to the electronic degree of freedom 
(DOF), the proton reaction coordinate x, and the proton donor-acceptor separation Rda, 

AG 

H s = hAa x + —a z + H a \a)(a\+H b \b)(b\ , (2) 

where \a) and \b) denote the electron donor and acceptor states, a x = \a)(b\ + \b)(a\, and 
a z = \b)(b\ — \a)(a\. H a and H b are defined as 

H ab = ^ + ^- + -m H co 2 H (x±-±n-) +-Mn 2 R 2 . (3) 

' 2m H 2M2 H \2 2J2 v; 

Here, p and x, P and R are the momenta and coordinates of the proton and proton donor- 
acceptor DOFs, respectively; m# and cu#, M and Q are the corresponding mass and frequen- 
cies; R is defined as the difference between the donor-acceptor separation and its equilibrium 
value R = Rda — R-da'i k * s ^ ne coupling coefficient between the proton equilibrium position 
and the proton donor- acceptor separation R, and ±(rf + kR)/2 is the equilibrium position 
of the proton coordinate on the donor and acceptor surfaces. 

The potential energy surface presented in Eq. (J3J) leads to shifted proton equilibrium 
positions on the donor and acceptor surfaces as a function of the proton donor-acceptor 
distance R. As the x and R motions are coupled differently on the \a) and |6) surfaces, 
there is a rotation of the corresponding normal modes (see also Fig. [T|), which is the DRE 



introduced in the previous Sec. HI In the case o 



the electron/proton donors and acceptors, 19M21 



a wide 



25 



y used linear model for the proton and 



45| k — 1. Without loss of generality, 



we will assume such case and drop k in the following derivations. 



We now briefly discuss the relevance of the DRE in PCET reactions. The V-shaped 
x-R potential energy surfaces for the donor and acceptor states as shown in Fig. 1 can 
actually be found in many previous publications studyi ng p roton transfer (PT) (e.g., Fig. 
1 in Ref . J46| ) and PCET reactions (e.g., Fig. 4 in Ref . [47] ) . as an effect that the proton 
equilibrium positions shift in different direction on the donor and acceptor surfaces as a 
function of the proton donor-acceptor distance. Under the harmonic approximation, these 
potential energy surfaces display a rotation of the normal modes that can be described by 
the DRE. The parameters for a DRE model can thus be obtained by analyzing the donor 
and acceptor potential energy surfaces. 

We also assume that the solvent DOFs couple only to the electronic DOF, and the bath 
Hamiltonian can be written as 

n-t^ + W'*) < (4) 

where Xj, Pj, rrij, Uj are the coordinates, momenta, masses, and frequencies of the harmonic 
bath modes. 

The bath DOFs are assumed to couple linearly with the electronic DOF, and the system- 
bath coupling term is given by 

N 

Hbs = ~J2 c i x i az ■ ( 5 ) 

3=1 

The system-bath coupling is usually characterized using the spectral density J(oj) defined 

as 

N c 2 

■^lE^e--*)- «■> 

3=1 J J 

In general, the proton DOF also couples to the environmental DOFs, which leads to a 



reorganization energy associated with proton motion, 12H14] as well as vibrational energy 



relaxation and dephasing. [48|, |49] Although we did not consider such coupling in this paper, 
it can certainly be incorporated by extending the above model Hamiltonian described in 
Eqs. flUED. 



B. Fermi's golden rule 



In the electronic nonadibatic limit where hA is small, first order perturbation can be 
applied, and the rate constant can be calculated using Fermi's Golden Rule: [28] 



fc = A' 



e-^C^dt 



where the correlation function C(t) is defined as 

1 



C(t) = -Tr 



e -/9flb e *tfot e -*ffit 



Here, Z is the partition function of the donor state, Z = Tre P H °, 

Hq = H a + H B + 2_^ ' 
and 



CjXj , 



H x = H b + H B -J2c. 



j x j ■ 



(7) 



(8) 



(9) 



(10) 



Since the proton and proton donor-acceptor DOFs are decoupled from the bath modes, 
C(t) can be calculated as 



Cit) = ^Tr 



e -P H a e iH a t e -jrH b t 



C B (t) . 



(11) 



C B (t) on the right hand side of the above Eq. ffTTl) arises from the bath contribution and 
is defined as 



c B (t) 



1 



Tre-^ 30 



Tr 



— PH B 0pjiH B0 t — \H B \t 



(12) 



where H Bom = H B ± ^ . CjXj. It can be calculated analytically using the spectral density 
J (a;), resulting in the following equation:^ 



C B (t) = expi-- / du— ^ [coth(f3hu / 2) (1- cos ut) + i sin out] I . (13) 

As in the many previous theories for the PCET reactions, we apply the following approx- 
imation for C B (t), which can usually be obtained by a high temperature approximation and 
short-time expansion: 

At 2 _• A x 

(14) 



i4i 



C B (t) « e W *» 

where A is the electronic reorganization energy 

AJ(ou 



A 



du- 
ct KLU 



(15) 



6 



The first term on the right hand side of the above Eq. (TTT1) . J-Tr 



p — PH a „jr Hat „— jr Hbt 



describes the time correlation function of a two-mode DRE model, and can be calculated as 






~P H a e % Hat e -% H l> t 



— I dxdRdx'dR'(xR\e^ Ha{t+inp) \x'R')(x'R'\e-^ Hbt \xR) . 

Ax 



(16) 



Since H a and Hi are Hamiltonians of harmonic oscillators, the matrix elements can be 



calculated analytically, 50] and the correlation function can then be obtained using Gaussian 



integrals. Petal 
Refs. 



29 



31 



34 



s of such calculation can be found in many previous publications such as 



391 ] . and will not be presented in this paper. 



We further note that the harmonic x-R model was employed in this study in order to 
obtain an analytical expression for the time correlation function in Eq. ( ITT]) . Including 
;he anharmonic effects numerically in some of the widely used linear models (e.g., Refs. 



2l|,U, 



4l|) would also be straightforward, as calculating Eq. ( 1T6"j) in such cases only need 



to solve a two-dimensional Schrodinger equation, which can be done routinely on modern 
computers. We note that a similar idea of solving a two-dimensional Schrodinger equation 
for the proton and proton accepter-donor motion was proposed in Ref. 20 1. 



C. Rate constant within the BO approximation for the proton motion 



Due to the small mass of proton or deuterium, a BO approximation is often applied to 
separate the motion of the proton coordinates x and the proton donor-acceptor separation 
R. For a fixed donor-acceptor separation R, the proton coordinate can be quantized, and 
the vibronic states can be calculated as 

,2 



V 



2mH 



+ V a (x,R) 



\<j) jta (x,R)) = E jia \(j) jja (x,R)) 



(17) 



where a = a,b denotes the electronic donor and acceptor states, and j denotes the vibrational 
states for the proton DOF. According to Eq. (ED, V ajb (x, R) = l/2m H uJ H (x±d/2±R/2) 2 , and 
Ej jCe and \<pj iOC (x, R)) can be obtained from the eigenstates of displaced harmonic-oscillators. 
In the case of small electronic coupling considered in the above subsection, neglecting the 
non-BO coupling terms between the vibronic states in Eq. f|T7|) wil 
reaction as nonadiabtic transitions between a group of vibronic states. 



'ormulate the PCET 



|22j By further 



assuming an initial equilibrium distribution on the donor vibrational states, the total rate 



constant can be calculated as 



"* / , -*/i,a"'/ii/ > 



(18) 



fJ,V 



where P^ a is the equilibrium population of vibrational state \(j) on the donor potential 
energy surface; and k^ u is the transition rate constant from state |/x, a) to \v,b). The rate 
constant k^ u are then calculated using the effective system Hamiltonian involving states 
|/i, a) and \u, b) 



14 



22] 






*>r\<f>v) (lA*, a) {v, b\ + c.c.) + [-AG/2 + + \/2)Hu H \ |/i, a) (//, a| 

P 2 1 

2M 2 



frA< 



[AG/2 + (z/ + l/2)^ H ]|z/,6)(z/,6| 



(19) 



Since fluctuations of the the donor-acceptor coordinate R can cause changes of the overlap 



22 



23 



251 1 there is no genera 

nn 



exact ana 



14 



mm 



yti - 



5l|. 



integral {<j>n\<j) v ) by orders of magnitudes, [15j-|20j, 

cal expression for k^ u , and additional approximations are often employed 6|,[8 

We will briefly present the results in applying two widely used approximations to the model 

system presented in the previous subsections II. A and B in the case of nonadiabatic ET 

reactions. 

I — l 

In the rate constant expression originally proposed by Kutnetsov and Ulstrup,|5l| the 
-R-mode is treated statically, and the rate constant is obtained as a thermally average over 
the classical Boltzmann distribution of R: 

{AG,*, + A) 2 






A< 



7T 



\k B T 



cxp 



dRP(R)\S^(R)\ 2 , 



(20) 



AXk B T 

where AG^ U = AG + [y — fi)hu B , S^ U (R) = (<fi^(R)\<f) u (R)) is the Franck-Condon overlap 
of the vibrational wave functions, and P(R) = a/ MVt 2 /27rk B Te~ Am2R2/2kBT '. This static 
approximation has recently been applied to PCET reactions in enzymes by Klinman and 
coworkers. 



19 



231 ] For cases where hVt > k B T } the above Eq. ( 12"U|) can be extended to 



include t 
-R-mode, 



re quan tum effect of the -R-mode by using a quantum mechanic distribution for the 



18 



52] 






A- 



7T 



Xk B T 



cxp 



(AGV + A)- 
A\k B T 



dRP qm (R)\S„ u (R)\ 



(21) 



where P qm (R) = ^1/271 (R 2 ) e - R2 ^ R2 \ and (R 2 ) = h/2MQcoth{f3hQ/2). This extended 
UK expression will be used in later calculations for the static -R-mode approximation. 



Another approximation is to expand S^ U (R) around the equilibrium position 
an exponential function, which is widely used in PT 17|, 



53 



54| and PCET 



20 



R = using 



22] theories: 



S^{R) « ^(0)e 



-atuvR 



(22) 



The problem is now equivalent to an extended spin-boson model with exponential 
coupling. [55] k^ can then be calculated as 20j, [22J 



h ■ I A <? 



<ftexp[c£,(C fl (0) + ^(t))]e- 1 ^ 4 -^-^' 



where 



C R (t) 



n 



2MQ 



[coth(/37ifi/2) cos Qt - i sin Qt] . 



(23) 



(24) 



In later studies, we will denote Eqs. ( 1T8|) . ( 1231) . and f l24l . as the exponential coupling 
approximation. 

III. RESULTS 

A. Rate constant and KIE 

We now apply the above theories to calculate PCET rate constants and KIEs. Fig. [I] 
shows the donor and acceptor potential energy surfaces from the Hamiltonian H a and Hi, 
defined in Eq. ([3]). We will assume uh = 3000 cm -1 throughout this paper, the other 
parameters used in Fig. [TJare M = 100 amu, £1 = 100 cm -1 , and d = 0.45 A. It can be seen 
that when R is smaller than its equilibrium value R = 0, the distance between the energy 
minima for the proton coordinate on the donor and acceptor surfaces becomes smaller, which 
will lead to enhanced proton tunnelling. A key feature of the DRE, which is the rotation 
between normal modes on the donor and acceptor potential energy surfaces, can be clearly 
seen in Fig. [U 

Fig. |2] shows the PCET rate constant as a function of the driving force AG using the DRE 
model, where the proton coordinate x is coupled to i?-modes with different parameters. A = 
30 kcal/mol, A = 100 cm -1 , d = 0.45 A, and T = 300 K were used in the calculation. It can 
be seen that coupling to the i?-mode can significantly increase the rate constants by orders of 
magnitude, especially when the donor-acceptor mass is small. We also note that the PCET 
rate constants keep increasing with more negative AG < —A. The reason is that the overall 



reorganization energy is very large after taking account into the contribution from the proton 
coordinate, and a real turnover of the PCET rate can only happen at very negative AG (see 



the inset o 
coworkers 



Fi g, pi) . This effect was also observed in recent studies by Hammes-Schiffer and 



41 



56| where the explanation is based on vibronic transitions to high vibrational 
states. Such high activation barrier for ET indicates that the proton transfer should happen 
mainly through quantum tunnelling. 

The KIE is an important character of reactions involving proton transfer, which is defined 
as the ratio of the rate constant for hydrogen to that for deuterium. The dependence of 



PCET rate constants and KI 
theoretical studies. 



14 



22 



s on various parameters has been widely discussed in previous 



24J More specifically, the effects of different donor-acceptor 



parameters and driving forces have been investigated by Hammes-Schiffer and coworkers 



recently. 4l|, |56[ The main purpose of the calculations below is, however, to quantitatively 
assess the applicability of the BO approximation for the proton motion: the static R-mode 
and the exponential coupling approximations will be evaluated using the analytically solvable 
model Hamiltonian presented in Sec. II. 

Fig. |3] plots the PCET reaction rate constants for proton and deuterium as a function 
of the donor-acceptor motion frequency. Comparisons with the rate constant expression 
using exponential coupling approximation (Eqs. [18j [231 and El]), and the extended UK 
expression within the static i?-mode approximation (Eqs. [T8l and |2T]) are also shown. Three 
different donor-acceptor masses M = 100 amu, 20 amu, and 7 amu were considered. The 
other parameters are A = 30 kcal/mol, AG = -5 kcal/mol, d = 0.45 A, and T = 300 K. It 
can be seen that the rate constant increases significantly with the decrease of the i?-mode 
mass, where larger proton donor-acceptor fluctuation leads to larger enhancement of the 
PCET rates. The extended UK rate expression Eq. (121]) agrees well with the exact FGR 
result except for the small mass and high frequency cases (M=20 amu and Q >500 cm -1 , 
M=7 amu and Q >400 cm^ 1 ), which indicates that the dynamical effect of the .R-mode 
becomes important only at high frequencies. We can also see that the exponential coupling 
approximation Eq. ( 1231) starts to overestimate the rate constants and becomes inaccurate 
for small mass and small donor-acceptor frequency (Q < 200 cm -1 for M=100 amu, Q < 400 
cm -1 for M=20 amu, Q < 700 cm -1 for M=7 amu), with a steep rise of the rate constants 
for low donor-acceptor frequency Q. The reason is that, the range of fluctuation for R 
is quite large in such cases, and the exponential approximation for the vibrational wave 

10 



function overlap in Eq. (j2"2l becomes invalid. 

Fig. H] plots the KIE as a function of the donor-acceptor frequency calculated with 
different rate expressions. It can be seen that the KIEs are small at low frequencies, which 
means that the enhancement of the PCET rate due to i?-mode fluctuations is more significant 
for deuterium. The extended UK method agrees well with the FGR result except for the case 
of high proton donor-acceptor frequencies (Q > 500 cm -1 for M= 20 amu, and Q > 400 
cm -1 for M= 7 amu.). The rate constant expression Eq. (123]) has problems in cases of 
small mass and frequency. This problem is more severe when the donor-acceptor mass M 
is small (e.g., the left parts of Fig. 4(b) and (c)) as the exponential approximation in Eq. 
(I2"2"j) becomes inaccurate when the fluctuation of R is large. We note that the failure of the 
exponential overlap approximation in calculating KIEs at low donor-acceptor frequencies 
has been discussed recently, |41|] while it is more quantitatively characterized in FigJUin this 
study by comparing to the exact FGR result. 

In summarizing this subsection, we have shown that the model Hamiltonian including the 
DRE captures the main features of PCET reactions presented in many previous theoretical 
studies, ja, uA, l22M24| The advantage of the current model is that, it can include the effect of 
the proton donor-acceptor motion without the BO approximation for the proton motion, as 
well as further approximations for the vibrational wave function overlap. Our quantitative 
evaluation of the static i?-mode and the exponential overlap approximations could also be 
helpful in developing more accurate PCET theory in more general cases. 

B. Vibrational coherence at small reorganization energy 

Recently, abnormal temperature dependence of KIE were observed in experimental stud- 



ies of several systems. |57l.l58| where the KIE increases as the temperature increases, implying 



a larger apparent activation free energy for proton transfer. Theoretical studies has pointed 



out that this could be a phenomena associated with low reorganization energies. 25|, |4l|, |42 ] 



An interesting finding in the literature is that, when the reorganization energy is small (i.e. 



in the low friction regime), coherent tunneling may play a role in the ET dynamics. 43|, |44] 
It is thus interesting to investigate whether coherent tunneling is relevant to the abnormal 
temperature dependence of KIE in the case of small reorganization energies. 

We plot in Fig. [5] the dependence of the PCET rate constant as a function of the driving 



11 



force AG for A = 3 kcal/mol. The other parameters are M=20 amu, Q = 300 cm -1 , and 
d = 0.45 A. The rate constants show oscillations as a function of AG, and this oscillatory 
behavior is more pronounced for proton than deuterium. This effect is readily explained in 



the picture of the BO approximation, 41 



42 



in Ref. 



56| with the most recent explanation presented 



56| . For example, using the extended UK type of expressions: as the vibrational 



energy levels for proton motion is 3000 cm -1 , or approximately 8.6 kcal/mol, the small 
reorganization energy A = 3 kcal/mol is not large enough to bridge the gaps between each 
pairs of vibronic transitions, and peaks for rate constants occur at approximately AG^ + A ~ 
0. As the deuterium case has a smaller vibrational spacing, this also explains the less 
oscillatory behavior for deuterium. Based on the same reasoning, Hammes-Schiffer and 



coworkers have a 
of small A.[4l|,|42| 



so explained the abnormal temperature dependence of the KIE in the case 



Here we provide a complementary view to the above explanation that was based on the 
BO approximation to the proton motion. In the FGR rate expression Eq. (0), the real part 
of the integrand Gag(^) = Re[e~ l ~ t C(t)} can be regarded as the (normalized) reaction flux 
correlation function whose integral gives the rate constant. Within first order approximation, 
the time dependent acceptor population at short times can also be estimated as 



j t P a (t) « k(t)P d (0) 



(25) 



where P a is the acceptor population, P d is the donor population, and 

k(t) = 2A 2 Re I dtC AG (t) . (26) 

Jo 

The normalized flux correlation function Gag(^) is plotted in Fig. [6] (a), for two different 
AG values at -3.5 kcal/mol and -6.5 kcal/mol. The corresponding population evolution is 
shown in Fig. [6] (b) by assuming Pd{0) = 1. We can see that due to the smaller damping 
effect associated with small A, the second and third peaks in the flux correlation function 
starts to contribute to the rate constant, and their contribution can be either constructive 
(for AG = -3.5 kcal/mol, the contribution from second peak is 40% from the first peak) or 
destructive, (for AG = -6.5 kcal/mol, the contribution from second peak is -70% from the 
first peak), which leads to peaks or valleys in the rate constant as shown in Fig. [5J We note 
that this just the vibrationally coherent tunneling effect observed in the previous numerical 



studies o: 



photo-induced electron transfer reaction involving a high frequency intramolecular 



mode,[44| although a slightly different model is used in the studies of Ref. [44 1 . 



12 



An interesting observation here is that, at AG values near the peaks and valleys of the rate 
constant curve in Fig. |5l the temperature dependence of the KIE shows different behavior. 
Fig. [7] plots the KIE as a function of temperature for AG = -3.5 kcal/mol (close to a peak), 
which decreases when the temperature increases; while for AG = -6.5 kcal/mol (close to a 
valley), the KIE shows abnormal behavior as it increases with increasing temperature. The 
dashed lines shows the results by including only the first peak contribution from Gag(^), 
which shows the normal temperature dependence in both case. Comparing the solid and 
dashed curves thus shows that the contribution from the oscillations in Gag(^) play an 
important role in the abnormal temperature dependence for AG = -6.5 kcal/mol. This 
observation also holds for AG values near other peaks and valleys. So we conclude that 
the abnormal temperature dependence of KIE in the case of small reorganization energy is 
related to coherent vibrational tunneling, which is a new explanation complimentary to the 
previous ones within the BO approximation of the proton motion. 

IV. CONCLUSIONS 

Using a model Hamiltonian based on the Duschinsky rotation effect, we were able to 
quantitatively investigate the applicability of the commonly used BO approximation for the 
proton motion in studying PCET reactions. The extended UK method using the static R- 
mode approximation was found to work reasonably well except for small donor-acceptor mass 
and high donor-acceptor frequencies. The exponential approximation for the vibrational 
function overlap was also tested in the model systems, and found to be inaccurate in cases 
of large proton-donor separation fluctuations when the proton donor-acceptor mass and 
frequency are small. This may suggest that more accurate analytical approximations are 
needed in such cases. 

In the case of small reorganization energies, we find that coherent vibrational tunneling 
can lead to oscillations of rate constants with respect to the driving force AG, and different 
temperature dependence of KIE at AG near the peaks and valleys on the rate-driving force 
curve in Fig. [5] Our explanation of these results based on coherent tunnelling is comple 



mentary to existing theories within the BO approximation for proton motion. 25|, |4l|, l42l. l56| 
As the vibrational coherent tunnelling is sensitive to couplings of the proton coordinate 
to the environmental DOFs, 44J extending the current study to incorporate such couplings 



13 



would certainly be interesting. The current work could also be extended to investigate 



60[ where non-equilibrium dynamics involving vibra- 



photo-induced PCET reactions. [59 
tional motion is important. 
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FIG. 1: Potential energy surfaces for the electronic donor (left) and acceptor (right) states. 
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FIG. 2: PCET rate constants as a function of AG. Curves from bottom to top: No coupling to 
the i?-mode; M = 100 amu, M = 20 amu, and M = 7 amu. The other parameters are: A = 30 
kcal/mol, A = 100 cm" 1 , Q = 300 cm" 1 , d = 0.45 A, and T = 300 K. 
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FIG. 3: Rate constants as a function of the proton donor-acceptor frequency O, for proton (H) and 
deuterium (D). (a) M = 100 amu, (b) M = 20 amu, (c) M = 7 amu. The other parameters are A 
= 30 kcal/mol, AG = -5 kcal/mol, A = 100 cm -1 , d = 0.45 A, and T = 300 K. 
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FIG. 4: KIE as a function of the proton donor-acceptor frequency Q. (a) M = 100 amu, (b) M 
20 amu, (c) M = 7 amu. The parameters are the same as in Fig. [3l 
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FIG. 5: PCET rate constants as a function of AG for small reorganization energy A = 3 kcal/mol. 
The other parameters are: A = 100 cm -1 , M = 20 amu, S7 = 300 cm -1 , d = 0.45 A, and T = 300 
K. 
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FIG. 6: (a) Normalized flux correlation function for AG = -3.5 kcal/mol (solid line) and -6.5 
kcal/mol (dashed line), (b) Population of the acceptor as a function of time. The other parameters 
are the same as in Fig. [5j 
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FIG. 7: KIE as a function of temperature for (a) AG = -3.5 kcal/mol, and (b) -6.5 kcal/mol. Solid 
lines are the exact FGR result, while dashed lines show the results by including only the first peak 
of the correlation function as shown in Fig. [5] (a). The other parameters are the same as in Fig. 
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